# data from: https://movementdisorders.onlinelibrary.wiley.com/doi/full/10.1002/mds.28052?saml_referrer

########### Coding for Parkinson Dataset #############

#!/bin/bash

#Parkinson dataset requires manifest file, also demultiplexes samples

#QIIME import files using manifest file
qiime tools import \
  --type "SampleData[SequencesWithQuality]" \
  --input-format SingleEndFastqManifestPhred33V2 \
  --input-path /mnt/datasets/project_2/parkinsons/parkinsons_manifest.txt \
  --output-path /data/parkinsons/demux_seqs.qza

# Convert to .qzv for visualization
qiime demux summarize \
  --i-data demux_seqs.qza \
  --o-visualization demux.qzv


#NO demultiplexing step required

#Denoising and trimming using DADA2
qiime dada2 denoise-single \
  --i-demultiplexed-seqs demux_seqs.qza \
  --p-trim-left 0 \
  --p-trunc-len 251 \            #251 because of the 251 b.p read length, and high confidence  at the last base called
  --o-representative-sequences rep-seqs.qza \
  --o-table table.qza \
  --o-denoising-stats stats.qza


#Filtering out the non-PD individuals to compare read depth rarefaction in qiime2 viewer
qiime feature-table filter-samples \
  --i-table table.qza \
  --m-metadata-file parkinsons_metadata.txt \
  --p-where "Disease='Control'" \
  --o-filtered-table filtered_table.qza



####### Visualizing DADA2 outputs, converting to .qzv files for viewing on QIIME2 web viewer ###########
qiime metadata tabulate \
  --m-input-file stats.qza \
  --o-visualization stats.qzv

qiime feature-table summarize \
  --i-table table.qza \
  --o-visualization table.qzv \
  --m-sample-metadata-file parkinsons_metadata.txt
  
qiime feature-table tabulate-seqs \
  --i-data rep-seqs.qza \
  --o-visualization rep-seqs.qzv




######### Generating files for taxonomic analysis #############
# Variable region V4 was amplified, thus can use the following Silva database: silva-138-99-515-806-nb-classifier.qza

# Taxonomic analysis if you don't have to train the classifier
qiime feature-classifier classify-sklearn \
  --i-classifier /mnt/datasets/classifiers/silva-138-99-515-806-nb-classifier.qza \
  --i-reads rep-seqs.qza \
  --o-classification taxonomy.qza
  
# Taxonomy barplots visualization
qiime taxa barplot \
  --i-table table.qza \
  --i-taxonomy taxonomy.qza \
  --m-metadata-file parkinsons_metadata.txt \
  --o-visualization taxa-bar-plots.qzv

#Filter out the chloroplast and mitochondria sequences
qiime taxa filter-table \
  --i-table table.qza \
  --i-taxonomy taxonomy.qza \
  --p-exclude mitochondria,chloroplast \
  --o-filtered-table table-no-mitochondria-no-chloroplast.qza

# Generate a tree for phylogenetic diversity analyses (rooted-tree.qza)
qiime phylogeny align-to-tree-mafft-fasttree \
  --i-sequences rep-seqs.qza \
  --o-alignment aligned-rep-seqs.qza \
  --o-masked-alignment masked-aligned-rep-seqs.qza \
  --o-tree unrooted-tree.qza \
  --o-rooted-tree rooted-tree.qza 




######## Exporting .qza files into an export folder #########

qiime tools export \
  --input-path table-no-mitochondria-no-chloroplast.qza \
  --output-path parkinsons_export 

#convert feature-table.biom to .txt file
biom convert -i feature-table.biom --to-tsv -o feature-table.txt

qiime tools export \
  --input-path rooted-tree.qza \
  --output-path parkinsons_export  

qiime tools export \
  --input-path taxonomy.qza \
  --output-path parkinsons_export  